#-------------------------------------------------------------------------------
# Name:        <1-x> eep
# Purpose:
#
# Author:      Gazizov
#
# Created:     27-11-2012
# Copyright:   (c) gazizov 2012
# Licence:     <your licence>
#-------------------------------------------------------------------------------

from pylab import *
from scipy.special import *

m_e = 5.1099892811E-4
m_p = 0.939
Egth = 2*m_e*(1+m_e/m_p)
print Egth


kB = 8.617332478E-5  # eV K^{-1}
T0 = 2.72556	     # T_{CMB} in K
kT = kB*T0
w0kT = 1/kT

nnrm = 411/(2*zeta(3,1))
nnrm2 = -nnrm*w0kT**2
print ('kT= %e, w0/kT= %e , nnrm= %e , nnrm2= %e ') % (kT, w0kT, nnrm, nnrm2)


def I(y):
    return nnrm2*log(1-exp(-w0kT*y))

yy = zeros(101,float)
fy = zeros(101,float)

for i in range(101):
    yy[i] = 10**(-10.+0.1*i)
    fy[i] = I(yy[i])

plot(yy,fy)

xscale('log')
yscale('log')

show()


##l = 11.9
##print 2*zeta(3,1)/pi**2*l**3

quit()

alpha = 7.2973525698E-3
r_e   = 2.81794E-13    # radius of electron in cm
ar2   = alpha*r_e**2   # in cm^2
print  ar2, ar2*1E30
quit()
Eg,avr1x = loadtxt('Data/Sigma_1_X.dat',unpack='True')
Eg *= 1E-3
avr1x *= ar2*1E30

for j in range(len(avr1x)):
    print (' %e %e ') % (Eg[j], avr1x[j])

##quit()

Egp,avr1xp = loadtxt('Data/Sig_1_X.dat',unpack='True')


xscale('log')
xlim(1E-3,20.4)
ylim(1E-3,10.)
yscale('log')
plot(Eg,avr1x)
plot(Egp,avr1xp,'ro')
grid(True)
show()